DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY  CA  93943-5101 


0.5.  RAVAL  PC- 

MONTEREY,  C  LI    - 


COMPUTER  SIMULATION  STUDIES 
OF  COPPER  ATOMS  IN  110  CHANNELS  OF  COPPER  CRYSTALS 


*  *  *  *  * 


RENE  W.  LEEDS 


COMPUTER  SIMULATION  STUDIES 
OF  COPPER  ATOMS  IN  110  CHANNELS  OF  COPPER  CRYSTALS 


by 
Rene  W.  Leeds 
Lieutenant,  United  States  Navy- 


Submitted  in  partial  fulfillment  of 
the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE 

IN 

PHYSICS 

United  States  Naval  Postgraduate  School 
Monterey,  California 
1964 


Ueh,  ft. 


r 


LfSRARY 
VS.  NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CALIFC 

COMPUTER  SIMULATION     STUDIES 
OF  COPPER  ATOMS  IN  110  CHANNELS  OF  COPPER  CRYSTALS 

by 
Rene  W.  Leeds 

This  work  is  accepted  as  fulfilling 

the  thesis  requirements  for  the  degree  of 

MASTER  OF  SCIENCE 

IN 

PHYSICS 

from  the 

United  States  Naval  Postgraduate  School 


ABSTRACT 

The  interaction  of  land  5  kev  copper  atom  (primaries)  impacting  the 
Oil  surface  of  copper  crystals  was  studied   using  an  "n-body  type" 
computer  simulation  program.  Each  primary  energy  was   studied  using 
both  the  Bohr  and  Gibson  II  copper  on  copper  potentials.     Primary  pene- 
tration and  lateral  drift  or  spread  were  investigated  as  a  function  of 
impact  point. 

The  n-body  program  has  greater  penetration  values  than  the  binary 

13 
model  used  by  Robinson  and  Oen      .     Drift  or  spread  (defined  as  wander 

if  energy  loss  rate  is  low)    occurs  in  the  impact  region  of  transition 

from  hard  interaction  (near  a  lattice  atom)  to  soft  interaction  (open  or 

channel  regions).     Penetration    and  spread  contours  are  presented. 

The  guidance  and  assistance  given  by  Associate  Professor  Don  E. 

Harrison,,  Jr.  of  the  U.   S.  Naval  Postgraduate  School  in  this  study  is 
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1.   INTRODUCTION. 

The  subject  of  crystal  bombardment  by  energetic  ions  or  atoms  may 
be  discussed  from  either  external  effects  or  internal  interactions, 
Neither  phenomenon  is  well  understood,  nor  have  adequate  theories  been 
devised  to  explain  these  interactions.    One  external  effect  is  described 
by  the  sputtering  ratio,  which  is  defined  as  the  ratio  of  ejected  atoms 
from  a  crystal  per  impact  atom.   A  more  general  term  for  crystal  modi- 
fication caused  by  incoming  energetic  atoms  is  radiation  damage. 

When  incoming  atoms  or  ions  (called  primaries)  penetrate  the  crystal, 
a  series  of  collisions  may  take  place  in  which  the  primary  energy  is 
distributed  to  the  crystal  atoms.    Often  crystal  atoms  leave    their  sites 
and  become  moving  atoms  called  secondaries.    The  orientation  of  the 
crystal,  relative  to  the  primary  beam  direction,  is  a  critical  parameter 
in  any  description  of  the  interactions  that  take  place.    Deeper  penetrations 
occur    when  primaries  enter  more  open  or  transparent  channels  in  crys- 
tals. 

The  primary  energy  effectively  categorizes  the  type  of  interaction 
which  will  occur.    Heavy  incoming  atoms  with  relatively  low  energies 
usually  interact  with  the  atoms  of  the  crystal  through  elastic  collisions. 
Such  atoms    deviate  markedly  from  straight  line  trajectories.    High  energy 
primaries  interact  by  ionization  and  electron  excitation  of  the  crystal 


and,  as  a  result,  lose  little  energy.    These  atoms  usually  follow 

1* 
relatively  straight  line  trajectories. 

This  study  is  concerned  with  the  ranges  and  spread  of  lower  ener- 
gy primaries  as  they  dissipate  energy  through  internal  elastic  inter- 
actions in  the  crystals.    Range  is  defined  as  the  depth  of  penetra- 
tion of  the  primary  from  the  surface.    For  purposes  of  this  study  total 
path  distance  and  vector  range  from  impact  point  to  a  specified  point 
are  not  needed.      Spread  is  defined  as  the  radial  displacement  of  the 
primary  from  the  axis  of  its  original  impact  line. 

These  introductory  sections  provide  a  background  in  the  recent 
experimental  and  machine  calculation  studies  of  low  kev  atoms  and 
ions  in  crystals.    The  topics  of  sputtering  and  basic  radiation  dam- 
age are  not  primary  objectives  of  this  thesis.    However,  several 

good  sources  of  information  are  recommended"    a  basic  text  for  intro- 

2 
ductory  study,  by  Dienes  and  Vineyard,     and  more  advance  articles 

by  Wehner,  et.  al.  ,     Nelson,  et.  al.  ,      '     and  Harrison,  et.  al.    . 

Although  recent  experimental  range  studies  have  been  performed 

with  heavy  ions  bombarding  single  crystal  targets;  the  nature  of  these 

experiments  and  the  lack  of  knowledge  concerning  potential  functions, 

make  the  comparison  between  laboratory  and  simulation  data  diffucult. 

Until  recently,  the  only  theoretical  potential  functions  available  were 

*    All  footnotes  refer  to  the  Bibliography 


those  from  homonuclear  atomic  interaction.    Thus,  the  computer  pro- 
grams have  used  primaries  and  crystals  of  the  same  type,  while  in 
laboratory  experiments  the  primaries  and  crystals  have  been  differ- 
ent substances. 

Virtual  leaders  in  the  field  of  experimental  range  studies  are  J.  A. 
Davies  and  his  co-workers  at  Chalk  River,  Ontario.    Their  technique 
is  both  sensitive  and  reliable,  but  requires  that  radioactive  gas  ions 
be  used  as  primaries.    The  approach  entails : 

.  .  .bombarding  an  optically  flat  aluminum  target  with  a  mono- 
energetic  beam  of  radioactive  ions,  and  then  measuring  the 
depth  of  penetration  of  the  incident  beam  by  dissolving  succes- 
ive  uniform  layers  of  aluminum  from  the  target  surface  and  measur- 
ing the  amount  of  radioactivity  in  each  layer.     The  ion  bombard- 
ments were  carried  out  using  an  electrostatic  accelerator  spec- 
ially designed  for  the  purpose.    The  technique.  .  .  (dissolving  ex- 
tremely thin  layers  of  known  thickness).  .  .consists  essentially 
of  two  steps  :electrochemical  oxidation  at  constant  voltage  in 
aqueous  ammonium  citrate,  followed  by  chemical  removal  of  the 
anodic  oxide  film  in  a  phosphoric  acid  chromium  trioxide  solution. 
Due  to  the  highly  protective  nature  of  the  anodic  oxide  film,  this 
process  permits  highly  reproducible  surface  layers  of  metal  as 
thin  as  37A  to  be  removed.     .  .  .The  process  may  be  repeated  as  of- 
ten as  desired,  so  that  it  is  possible  to  obtain  the  complete  range- 
distribution  curve  from  a  single  bombardment. 

Some  of  the  earlier  polycrystaline  range  studies  performed  by  Davies, 

et.  al.  included  .7  to  60  kev  Na24  ions  in  aluminum,     and  2  to  600  kev 

or  8 

Kr°°  ions  in  aluminum  and  tungsten.      All  of  the  range  distribution  curves 

exhibited  "tails"  of  deep  penetration.    This  phenomenon  was  attributed  to 

channeling  of  ions  in  preferred  directions  in  the  crystal.     Davies,  et.  al. 


stated, 

The  magnitude.  .  .  (of  crystal  lattice  effects).  .  .particularly 
emphasizes  the  need  for  using  either  single  crystals  or 
amorphous  materials  in  future  range  experiments .  ^ 

q 
In  a  more  recent  study  the  same  group     studied  single  crystals  of 

tungsten  using  1  to  20  kevXe125    ions.    Tungsten  has  a  body  center- 
ed cubic  (b.c.c.)  structure  in  contrast  to  aluminum  and  Cu  which  have 
face  centered  cubic  (f.c.c.)  crystals.    General  results  of  the  single 
crystal  studies  were. 

. .  .crystallographic  effects  are  even  larger  than  those  found 
previously  in  Al;  furthermore  the  results  are  in  qualitative 
agreement  with  the  theoretical  prediction  for  b.c.c.  lattice, 
i.e.  the  (100)  and  (111)  are  the  most  favoured  directions  for 
iteanneling  and  the  (110)  and  (112)  are  less  favoured.    This 
contrasts  with  the  f.c.c.  structure,  where  both  theory  and 
experiment  find  (110)  the  most  favoured  and  (111)  one  of  the 
less  favoured.9 

Lutz  and  Sizmann       bombarded  copper  (f.c.c)  single  crystals  with 

Kr  ions  of  10  to  140  kev.  energy.    Their     results  agreed  in  general  with 

Davies1  single  crystal  data,  since  they  showed  that  110    channel  shots 

had  deepest  penetration.    Further,  the  range  distribution  curves      showed 

o 

the  characteristic  tail  of  deep  penetration  (up  to  4000A)  for  110,  111, 

and  100  channels.    Their  25  kev  and  140  kev  range  distribution  curves 
(as  well  as  some  of  Davies'  curves)  are  shown  in  figure  1.    The  experi- 
mental procedure  used  by  Lutz,  et.  al. 10    employed  the  removal  of  thin 
layers  of  bombarded  crystal  by  controlled  sputtering  with  non-radio- 
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active  krypton  ions.    The  radioactivity  was  then  measured. 

In  conjunction  with  the  monocrystal  experiments  described  above 
and  computer  simulation  studies  to  be  discussed,  Lehmann  and  Leib- 
fried  have  studied  the  behavior  of  primaries  in  channels  analytical- 
ly. Their  work  used  copper  primaries  in  copper,  as  have  the  majori- 
ty of  computer  programs.    Their  results,  oriented  toward  long  range 

penetrations,  in  part  were: 

Two  potentials  are  used:  an  exponentially  screened  Coulomb 
potential  after  Bohr,  used  also  in  the  machine  calculation  and 
thought  to  give   an  adequate  description  for  relatively  high  en- 
ergies and  small  interatomic  distances;  and  a  purely  exponen- 
tial potential  after  Born-Mayer,  better  suited  for  relatively  low 
energies  and  large  atomic  distances.    The  maximum  ranges  are 
very  large,  for  10  kev  in  the  order  of  103  lattice  parameters  for 
the  Born-Mayer  potential  and  up  to  10'  for  the  Bohr  potential. 
Presumably,  the  Born- Mayer  potential  is  a  better  description  for 
these  events.     The  investigation  is  confined  to  motions  near 
the  channel  axis . 


2.     COMPUTER  SIMULATION  STUDIES 

This  project  studies  the  interactions  of  energetic  atoms  in  crys- 
tals through  the  use  of  a  digital  computer  simulation  program .  This 
section  will  present  a  general  orientation  of  the  use  of  computer  sim- 
ulation programs  and  will  discuss  the  current  work  in  this  area.    (A 
more  detailed  discussion  of  the  program  and  model  used  in  this  study 
will  be  presented  in  section  four  and  in  Appendix  II.) 

Simulation  programs  for  any  system  have  the  following  three  es- 
sential ingredients.     (1)  The  system  or  phenomenon  has  a   status. 
This  status  is  described  in  terms  of  attributes.     In  the  case  of  a  crys- 
tal, such  attributes  as  position,  velocity,  mass  and  energy  are  essen- 
tial.    (2)  The  status  is  changed  or  modified  by  an  event.     Events  cause 
the  atoms  to  leave  their  sites,  gain  or  lose  energy,  etc.     (3)  The 
third  essential  element  of  a  simulation  program  is  time.    The  status  and 
events  are  strongly  intertwined  in  a  time  continuum.     The  control  of  the 
occurance  of  events  is  critical. 

Three  basic  approaches  to  simulation  programs  of  atom-crystal  in- 
teractions that  have  been  used  to  date  are:  (l)random  target  atoms  (i.e. 
lattice  effects  are  neglected),     (2)  binary  collision  in  a  lattice,  and 
(3)  many  body  collision  in  a  lattice.    The  binary  collision  assumes  in- 
teraction with  one  atom  at  a  time,  while  the  many  body  approach  (called 

n-body)  simultaneously  accounts  for  forces  from  several  crystal  atoms. 

12 
One  of  the  earliest  computer  studies  was  made  by  Gibson,  et.  al. 
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Their  model  placed  the  copper  atoms  on  sites  in  the  f.c.c.  crystal 
and  used  an  n-body  approach.     To  simulate  an  infinite  crystal,  add- 
itional forces  were  added  to  the  surface  atoms.     Because  the  poten- 
tial functions  and  their  corresponding  force  functions  were  not 
known  precisely,  Gibson,  et.al.   used  "a  simple  central  difference 
procedure.  .  .which  gives  reasonable  accuracy.  "*■*    Their  techniques 
used  f=ma  in  an     iteration    procedure.    This  study  showed  evidence 
of  chains   in  the  100  and  110  directions.  A  chain  is  defined  as  the 
propagation  of  energy  along  close  packed  rows  of  atoms  in  a  crys- 
tal without  mass  transport.     Interatomic  potential,   surface  effects, 
lattice  defects,  and  other  low  energy  (<400  ev)    phenomena  were 
studied. 

The  random  lattice  mentioned  above  was  used  in  an  earlier   study 
by  Oen,  et.  al.     for  copper  primaries  in  solid  copper.     Their  pro- 
gram assumed  independent  binary  collision  (accomplished  by  an  im- 
pact parameter  restriction  of  one  half  the  nearest  neighbors  distance), 
classical  scattering,  and  the  Bohr  hard  sphere  approximation  to  the 
Bohr  potential.     They  concluded: 

It  is  found  that  neither  the  hard  sphere  approximation  nor  the 
inverse  r  squared  approximation  to  the  Bohr  potential  is  partic- 
ularly good.     To  obtain  correspondence  with  experimental  re- 
sults it  is  found  that  the  Bohr  screening  length  must  be  in- 
creased as  the  atomic  number  of  the  interacting  atoms  increas- 

PQ       1 


A  more  recent  program  by  Robinson  and  Oen       (hereafter  referred 
to  as  RO)  used  the  binary  approach  again,  but  replaced  the  random 
distribution  with  a  copper  f.cc.  lattice.    This  program  studied  sev- 
eral potentials;    the  Bohr,  eroded  Born- Mayer,  and  truncated  Born- 
Mayer»     The  Born-Mayer  potential  selected  by  RO  had  the  same  para- 
meters as  the  Gibson  II  potential,  but  incorporated  four  different 
truncation  values. 

Significant  among  their  conclusions  were  the  following: 

(1)  Ranges  of  primaries   were  strongly  dependent  on  direction, 
and  in  f.cc.  crystals  the  order  of  range  was  (01 1)>  (001)>  (111)?« 
isotropic. 

(2)  Larger  range  values  were  due  to  many  glancing  collisions  that 
confined  the  primaries  to  regions  of  low  potential.     This  phenomenon  of 
deep  penetration  is  called  channeling. 

(3)  The  combination  of  Bohr  potential  and  1000  ev  primaries  im- 
pacting in  the  (Oil)  direction  produced  several  primaries  that  moved 
from  channel  to  channel.    That  is,  RO  found  that  the  primaries  showed 
large  spread  as  they  penetrated  the  crystal.     Such  spread  was  attributed 
to  the    nature  of  the  Bohr  potential,  feee  figure  2,  which  is  a  projection 
of  their  (Oil)  channel  trajectories  onto  the  (Oil)  plane.     Note  that  three 
of  the  four  tracks  show  large  amounts  of  spread.) 
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A  more  recent  n-body  model  was  used  by  Gay  and  Harrison 
(hereafter  referred  to  as  GH).     They  used  an     iteration  procedure, 
but  one  that  has  important  differences  from  Gibson's  approach 
(see  Section  4).    Further,  the  crystal  size  was  smaller  than  Gibson's, 
and  the  program  was  simplified  to  reduce  computer  time.    A  compari- 
son between  binary  collision  characteristics  and  n-body  collision 
approximations  was  performed.     In  general  terms,  the  report  showed 
that  binary  collision  approximations  approach  true  n-body  collisions 
for  small  impact  parameters.     But  at  larger  impact  parameters  the 
binary  collisions  have  less  energy  transfer  from  primary  to  target,  re- 
duce the  primary  scattering  angle  and  slightly  increase  the  recoil 
angle.     Their  work  provides  specific  quantitative  results  for  interac- 
tions in  terms  of  energy  exchanged,  scattering  angles,  and  recoil 
angles. 

A  modified  form  of  the  GH  model  was  used  in  this  study. 


3 .     STUDY  OBJECTIVES 

This  study  of  energetic  atoms  in  crystals  was  performed  using  an 
n-body  model  approach,  incorporated  in  a  computer  simulation  pro- 
gram.    Specifically  the  study  objectives  were: 

(1)  The  verification  and  more  explicit  delineation  of  the  spread 

13 
phenomenon  found  by  RO,        (figure  2).    As  mentioned  earlier,  three 

of  the  four  copper  primaries  shown  in  the  figure  drifted  considerably 

in  the  lateral  direction  as  they  penetrated  the  crystal. 

(2)  The  determination  of  the  range  distribution  curves  for  copper 
atoms  penetrating  the  (Oil)  surface  of  f.c.c.  copper  crystals.     Spec- 
ific attention  was  placed  on  the  effect  of  impact  point  on  primary 
range  and  spread. 
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4o     SIMULATION  MODEL 
A.     Model  Discussion 

The  original  purpose  of  the  GH  program  was  a  comparative  study  of 
binary  versus  n-body  collision  processes .    Their  model  used  a  63  atom 
f.cc.  microcrystalliteo    These  atoms  were  capable  of  movement  but 
could  not  dissipate  energy  except  by  further  collision     No  surface  ef- 
fects were  allowed  ;  which  was  clearly  justified  by  the  short  time  span 
of  the  program  .    An  outer  shell  of  immovable  atoms  was  used  to  simu- 
late a  continuation  of  the  lattice  and  test  for  complete  event  contain- 
ment.    (This  shell  is  eliminated  in  the  current  program), 

The  objectives  of  the  GH  study  were  oriented  toward  analysis  of  the 
primary  atom  and  one  "target"  atom.    Concentration  on  the  primary's 
history  allowed  a  streamlining  of  the  program  that  removed  unneeded 
operations  and  reduced  program  run  time,    Essentially  the  computation 
and  recording  of  the  motion  of  the  remaining  crystal  atoms  were  deleted. 
This  deletion  was  possible  since  the  path  of  a  primary  is  rarely  exposed 
to  feedback  from  target  interactions.     In  analogy,  if  the  path  of  radia- 
tion damage  is  represented  by  an  expanding  "tree  of  interactions",  then 
the  primary  is  usually  spearheading  the  top  branches  „ 

Another  result  of  the  concentrated  study  on  a  primary  and  a  single 
target  was  the  reduction  in  crystal  size.     It  must  be  made  clear  that  the 
reduction  in  crystal  size  has  been  carefully  studied  in  earlier  work  by 
Harrison.         Considering  the  uncertainties  of  the  crystal  potentials  and 
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the  use  of  iteration  techniques,  the  neglect  of  more  distant  crystal 
atoms  is  not  significant  in  regard  to  primary  performance.     The  crystal 
reduction  was  carried  out  experimentally  through  trial  runs  on  the  com- 
puter.    The  size  was  reduced  until  the  program  run  time  was  optimized, 
but  iteraction  results  still  correlated  closely. 

The  unit  crystal  used  in  this  study  regenerates  itself  in  one  of  two 
ways.     The  first  form  constructs  the  crystal  in  the  y  direction  only 
(figure  3)  while  the  second,  slower  program  can  construct  the  crystal  in 
any  direction  (figure  4).    As  a  result  of  these  features,  the  small  unit 
crystal  is  repeatedly  remade  in  front  of  the  primary,,     The  fast    running 
program  is  used  in  those  cases  where  the  primary  is  not  expected  to  ex- 
perience severe  collisions  nor  deviate  greatly  from  its  original  track. 
This  type  of  program  is  called  Chan2.     The  second  program  is  used  in 
areas  when  the  primary  does  experience  large  amounts  of  spread  and  is 
designated  Chan  3.     Figures  3  and  4  indicate  the  unit  crystal  sizet  and 
their  orientation.    All  primaries  hit  the  Oil  surface  orthogonally  (i.e.  in 
the  y  direction).     The  crystal  is  positioned  by  the  program  with  Oil  sur- 
face in  the  x-z  plane  and  with  the  positive  y  direction  coinciding  with 
the  primary  track  at  impact.    A  majority  of  the  runs  were  made  in  a  tri- 
angular impact  area.     Several  check  runs  were  made  in  the  adjecent  triangle 
and  a  few  elsewhere.     Detailed  discussion  of  the  impact  triangles  is  given 
in  the  next  section. 

o 

The  basic  building  block  of  a  copper  crystal  is  a  cube  measuring  3.614A 

on  a  side.    The    simulation  program  uses  units  of  length  of  1.807A  called 
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lattice  units . 

The  program  is  written  in  Fortran  and  Symbolic  Fortran  (an  assembly 
language  closely  related  to  machine  language).    All  runs  were  conducted 
on  a  Control  Data  Corporation  1604  computer  which  has  a  storage  capacity 
of  over  32,000  words  and  six  index  registers.    A  CDC  160  computer  was 
used  to  automatically  plot  graphic  projections  of  the  primary  trajectories. 
Refer  to  figure  5  for  an  example  of  such  tracks  projected  onto  the  011  plane. 

The  program  constructs  the  lattice  (i.e.  assigns  positions  to  each  atom), 
commences  the  run,  at  designated  time  intervals  computes  the  forces  on  the 
primary  and  secondaries,  adjusts  their  motion  accordingly,  prints  required 
information  at  set  intervals,  and  "shuts  down"  when  any  of  three  conditions 
exist.    The  program  terminates  when  (1)  the  energy  of  the  primary  goes  be- 
low 25ev,  (2)  a  preset  time  limit  is  exceeded  (measured  in  program  cycles; 
10  cycles  « time  for  the  primary  to  travel  one  lattice  unit),  or  (3)  the  primary 
leaves  the  lattice  (necessary  for  Chan  2  programs  only).    The  determination 
of  this  cycle  factor  received  considerable  study  by  Gay   "  and  was  optimized 
experimentally. 

As  may  have  been  evident  in  the  brief  program  discription  above,  the  crit- 
ical phas-3  of  the  program  is  the  method  of  computation  of  the  primary  path  as 
it  experiences  crystal  forces.     In  general,  the  "events"  of  this  simulation 
program  are  the  primary  crystal  interactions. 

The  interaction  routine  used  by  GH  is  still  used  in  this  model.     It  consists 
of  a  double  iteration  procedure  that  calculates  the  average  force  on  an  atom 
over  a  small  distance  of  travel.    This  force  is  used  to  advance  the  atom. 

m    » 
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Harrison       states; 

The  unbalanced  force  .  .  „  (on  the  primary).  .  ,  is  an  average  force 
calculated  by  a  double  iteration  procedure  as  follows    (1)  assume 
an  atom  at  position  1  with  velocity  1;  (2)  calculate  the  total  force 
on  the  atom  as  a  result  of  all  the  other  atoms  in  the  lattice  (this 
means  normally  only  about  8-10  nearest  atoms.  .  .);  (3)  call  this 
calculated  force,  force  1,  and  use  the  equation  of  motion  to  move 
the  atom  to  a  temporary  position,  position  2;  (4)  now  repeat  the 
force  calculations  for  position  2,  call  this  force  2;  (5)  go  back  to 
positional,,  and  use  the  average  of  force  1  and  force  2  to  move 
the  atom  to  a  new  position,  position  3.     Procedures  1  through  5 
constitute  one  "time  step".-'-4 

The  equation  of  motion  used  in  the  iteration  procedure  differed  in  two 

12 
ways  from  the  equation  used  by  Gibson       et.  al.    Although  derived  in  a 

different  manner,  the  GH  equation  can  be  compared  more  easily  to  the 

Gibson  form  by  starting  with  the  latter  and  indicating  the  modifications. 

The  Gibson  form  of  the  equation  (in  the  x  direction  only)  is; 

Xi  (t+  At)  =Xi  (t)  +  At[vi  (t+  At/2)  +  m"1.  F  (t)  .  A  t j  , 
where  X^^  =  initial  X  coordinate,  t  =  time,  At  =  time  increment  for  change 
of  position,  vi  =  initial  velocity,  m  =  mass,  and  f=  the  initial  force. 
Because  At  is  a  small  number  compared  to  t  the  term  Vi  (t  +  At/2)  may  be 
expanded  in  a  Maclaurin  series,     Neglecting  all  terms  after  the  first  ,  the 
following  form  is  obtained-  Xj  (t  +  At )  »  xi  +  At  [  Vi  (t)  +  At/2.  F(2)/Mj 
In  the  GH  derivation ,  the  force  term  is  an  average  force  where 

:        Fi+  Ff 
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and  the    following  expression  results t     -       -> 
Xi  (t  +  At  )  =  Xi  +  A  t  [Vi  (t)  +  41*  TJ^J 

Thus,  the  GH  expression  resembles  the  Gibson  expression  but  utilizes 
a  more  realistic  force  term„ 
B.     Potential  Functions 

Two  factors  dominate  all  other  considerations  in  the  simulation  of  pri- 
mary interaction  in  crystals .    The  first  is  the  general  type  of  model  attempt- 
ed. Specifically,  does  the  program  use  binary  or  n-body  collisions?     Is  a 
crystal  or  a  random  atom  arrangement  used?    The  second  factor  is  the  se- 
I:    tion  of  a  potential  function.     Many  functions  have  been  proposed;  none 
have  been  singled  out  as  best.    Until  good  correlation  between  experiment 
and  computer  simulation  is  made,    the  proper  potential  function  will  remain 
uncertain.     In  a  sense  this  report  and  all  such  simulation  studies  are  an  at- 
tempt to  locate  the  correct  function  for  a  particular  substance. 

Two  basic  types  of  potential  function  were  used    in  this  study"    Born- 
Mayer  and  Bohr.    The  Born-Mayer  function  used  has  three  parameter  varia- 

12 
tions  introduced  by  Gibson  et.  al.        and  shown  below: 

0=Aexp[-p  (r-r0 )  /rj    , 
where  r  =  distance  in  lattice  units, 

r    =  nearest  neighbors  distance  (2.555A    for  copper), 

o 

10=  lattice  unit  (1.804A    for  copper), 
p   =  input  parameter,  and 
A    =  input  parameter  (ev) . 
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The  form  shown  above  is  convenient  for  computer  calculation.    The 
three  sets  of  parameters  selected  by  Gibson  are: 

Potential                                            A(ev)  p 

I                                                  .0392  16„97 

II                                                  .051  13.00 

III                                                   .1004  10.34 
To  speed  calculation,  the  Born- Mayer  force  functions  were  truncated  at  f=10 

o 

newtons  (corresponds  to  approximately  3A). 

In  similar  format  the  Bohr  potential  function  is: 

0  b  =L1/rJ   e  x  p  {Sn  cb  ab  -  r/ab  ), 
where  r  =  distance  in  lattice  units, 

ab=  input  parameter  in  lattice  units,  and 

cb=  input  parameter  in  ev. 

4 
The  two  parameters  are  i  cb  =  9.94x10  ev,  and  ab=  .06741u.    These  para- 
meters are  the  ones  used  in  the  RO  report.    Further,  the  corresponding  force 
truncation  value  was  used  so  that  the  Bohr  potential  functions  and  force  func- 
tions in  this  study  are  identical  to  Robinson's .   (The  truncation  of  the  force 
function  occurs  at  f  =  3.9x10       newtons). 

The  lattice  has  zero  kinetic  and  potential  energy  prior  to  primary  impact. 
In  other  words  no  vibrational  energy  is  present. 
C.     Binary    vs  n-body  model 

This  study  made  use  of  an  n-body  simulation  program  in  contrast  to  the 

1    13    17 
simplified  binary  collision  approach  of  Robinson,      '      '       et„  al.    The  in- 
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herent  danger  of  superimposing  independent  collisions  is  that  simultaneous 
collisions  will  not  be  properly  accounted  for.    For  purposes  of  fast  comput- 
er runs,  the  binary  program  has  the  greatest  advantage,  but  its  justifica- 
tion requires  better  experimental  verification  than  is  yet  available.    The 
n-body  program  is  inherently  more  accurate  because  of  its  more  thorough  in- 
teraction routine,,    Given  proper  parameters,  it  should  yield  the  more  accurate 
data. 
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5  o     EXPERIMENT  PROCEDURE 

Computer  simulation  studies  are  experiments  in  the  accepted  sense  of 
the  word.    They  are  trials  made  to  confirm  or  disprove  some  suggested 
phenomenon,  and  are  subject  to  most  of  the  usual  experimental  inaccuracies . 
Systematic  errors  occur  because  of  improper  model  simulation  or  necessary 
oversimplifications  of  the  model.    Random  errors  caused  by  computer  mal- 
function should  be  rare„     Provided  the  program  and  its  input  data  are  not 
modified,  the  computer  should  reproduce  results  perfectly.     (Refer  to  subsec- 
tion C . ) . 

Simulation  studies  are  sensitive  to  parameter  variation  to  the  extent  that 
some  computations  can  become  meaningless  if  not  properly  controlled.    For 
example,  the  combination  of  high  energy  primaries,  Bohr  potential,  and  sel- 
ected impact  points  produce  an  interaction  so  weak  that  the  program  cannot 
maintain  the  energy  check.    As  a  result  the  primary  could  penetrate  infinite- 
ly.   Such  extreme  situations  have  been  avoided  and  are  not  included  in  the 
data  of  the  report,  but  they  do  point  out  the  danger  of  blind  faith  in  computer 
"answers" . 

The  remainder  of  this  section  will  present  the  experimental  approach  used. 
A.    Determination  of  impact  area  and  sample  size 

Symmetry  considerations  reveal  that  the  impact  rectangle  shown  in  figure 
3  or  4  represents  all  possible  combinations  of  impact  possibilities.  That  is, 
the  rectangle  is  a  "unit  area"  that,  if  repeated,  will  completely  describe  the 
Oil  surface.    Further,  with  one  exception,  the  rectangle  can  be  split  down 
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the  diagonal  into  two  triangles  that  have  point  correspondence  in  regard  to 
primary  interaction.    Consider  impacts  at  points  A"  and  C'(the  90  degree  an- 
gles).   As  indicated  in  figure  4,  atom  row  AB  is  nearer  the  surface  than  row 
DC  and  therefore  "sees"  the  primary  sooner.    Now,  a  primary,  or  "bullet", 
striking  the  crystal  at  A'  "sees"  these  top  plane  or  near  atoms  (row  AB)  at 

o 

a  distance  of  1.28A  .    However,  a  bullet  hitting  at  point  C'  sees  the  top 

o 

plane  atoms  (row  AB)  at  a  distance  of  1.807A  .    Obviously,  this  difference 
effects  the  first  interaction  and  subsequently  alters  the  bullet's  path  and 
future  collisions.     It  will  now  be  shown  that  this  difference  in  the  two 
triangles,  will  cause  a  mild  perturbation  in  range  and  spread  data,  and  is 
not  considered  serious  for  purposes  of  this  study „     (The  triangular  impact 
area  is  desirable  since  it  reduces  the  sample  size  necessary  by  one  half), 

A  graphical  correlation  analysis  of  two  critical  parameters  (range  and 
spread)  was  performed  to  insure  the  validity  of  the  impact  triangle.     In 
each  triangle  14  points  were  chosen  and  the  range  and  maximum  spread 
values  recorded.    All  points  were  obtained  using  lOOOev  primaries  and  a 
Gibson  II  potential.    Since  this  combination  has  a  strong  primary-crystal 
interaction,  discrepancies  between  the  two  triangles  should  be  maximum 
(figure  7).     In  general,  good  correlation  was  revealed.     Note  that  the  range 
correlation  values  show  less  correlation  for  deeper  penetrations.    This  is 
expected  since  those  primaries  that  impact  further  from  the  "corner  atoms" 
have  smaller  energy  loss  rates  and  more  chance  to  deviate.    The  spread 
terms  have  less  correlation  at  small  values.     However,  these  terms  are 
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sufficiently  small  so  that  .2  or  .3  lattice  unit  variation  between  triangles 
does  not  effect  the  spread  contour  diagrams  (to  be  discussed  later).    As  a 
result  of  these  considerations,  triangle  ABD  was  selected  as  the  primary- 
impact  area.    The  error  introduced  by  this  assumption  is  small  and  is  not 
known  precisely.     In  future  studies  one  area  of  further  investigation  should 
undoubtedly  be  a  more  thorough  check  of  the  rectangular  impact  area  and  its 
effect  on  results. 

From  the  grid    shown  in  figure  8,  14  impact  points  were  selected  and 
are  labeled  one  through  14  «>    For  some  of  the  runs  these  original  14  points 
were  supplemented    by  12  more,,    These  points  are  labeled  in  the  figure  as 
1A  through  12A.    The  supplemental  12  point  array  was  chosen  to  verify  the 
selection  of  the  basic  14  point  sample  size  .    Both  the  Gibson  11=1000  and 
5000ev  studies  were  initially  made  on  14  point  grids  and  later  supplemented 
by  12  more  data  points.    The  range  contour  surfaces  remained  regular  when 
the  new  points  were  added,  and  showed  no  anomalies  (figure  9),    Further, 
the  range  distribution  curves  for  the  14  point  and  26  point  grids  correlated 
very  closely  (figure  6),  and  established  the  soundness  of  14  point  grids. 
These  range  curves  were  constructed  on  an  area  rather  than  a  point  basis. 
That  is,  cross  sectional  areas  of  the  contour  surfaces  (figures  9  through  12) 
were  measured  at  specified  depths  of  penetration „    The  ratio  of  this  area 
to  the  total  triangle  area  gives  the  percent  of  primaries  stopped.     In  figure 
9  the  shaded  area  represents  those  primaries  stopped  within  100  lattice  units. 
The  ratio  of  this  area  to  the  total  area  provides  the  percent  stopped  at  100 
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lattice  units.     (The  range  curves  were  constructed  from  larger  contour 
diagrams  than  those  shown).    The  area  method  yields  more  accurate  pene- 
tration curves  than  would  a  point  system.    Overall  justification  of  a  14 
point  sample  is  simply  the  general  regularity  of  the  results. 
B.     Impact  points  A, B,C  and  Do 

In  addition  to  the  26  point  grid,  four  other  points  were  chosen  on  the 
triangle.     These  points  coincide  with  the  impact  points  in  figure  3  of  the 
RO  study  (three  of  the  four  points  showed  spread).     These  points  are 
shown  in  figure  8  of  this  report  and  are  labeled  as  A,B,C  and  D. 

As  an  initial  step  in  this  study,    the  parameters  of  the  RO  "A,B,C 
and  D  shots"  were  duplicated  as  closely  as  possible.     Later,  several 
other  combinations  of  bullet  energy  and  potential  were  used.    Table  1 
shows  the  number  of  grid  points  investigated  under  each  combination  of 
parameters.    As  will  be  shown  in  the  next  subsection,  many  of  these 
points  were  rechecked  with  additional  runs. 
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TABLE  1 


Impact  Points  Studied  By  Parameters 

RO13  Points 
A,B,C,D 

14  Point 
Grid 

12  Point 
Grid 

Symmetrical 

Triangle 

Potential 

Energy 
1000     5000 

Energy 
1000     5000 

Energy 
1000    5000 

Energy 
1000 

Bohr 

4 

Point  C 

14 

13 

8 

6 

- 

Gib  II 

4 

Point  C 

14 

14 

12 

12 

14 

Gib  III 

4 

- 

14 

- 

- 

- 

- 

C.    Reproduceability  of  runs. 

As  stated  earlier  reproduceability  problems  should  not  occur  in  com- 
puter studies.     Unfortunately,  the  program  -  computer  combination  used 
in  this  study  does  fail  to  reproduce  periodically.    As  a  result  a  large 
number  of  reruns  were  necessary.    The  sample  was  random,  but  was 
later  checked  to  insure  that  most  combinations  of  energy,  potential  and 
run  length  were  represented.    The  sample  size  was  approximately  15  per- 
cent.    It  was  found  that  the  probability  that  a  run  would  be  reproduced 
exactly  was  85  percent.    Further,  once  a  run  was  confirmed  by  a  rerun, 
the  probability  of  its  correctness  could  be  assumed  to  be  one.    The  lack  of 
perfect  reproduceability  has  not  been  determined.     However,  it  was  found 
that  those  runs  that  deviated  seemed  to  get  out  of  phase,  but  did  not  vary 
greatly  in  amplitude  or  distance  traveled.    Generally,  the  variation  of  am- 
plitude (i.e.   spread)  and  depth  of  penetration  were  less  than  a  few  per- 
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cent.     In  one  case  the  depth  of  penetration  varied  by  about  50  percent, 
but  the  depth  of  penetration  in  this  case  was  very  short  (  <8Lu)» 
D.     Program  outputs 

The  output  will  be  briefly  discussed  here,  so  that  subsequent  dis- 
cussions of  data  will  have  more  meaning.    The  program  as  used  in  this 
study  was  designed  specifically  to  examine  only  the  primaries'  per- 
formance.   Output  is  made  every  ten  cycles.    Ten  cycles  is  the  time 
necessary  for  the  bullet  to  travel  approximately  one  lattice  unit  (1.807A). 
Each  output  yields  the  following  data  about  the  primary- 

(1)  x,y,  and  z  distances  (measured  in  z  lattice  units)  from  the  im- 
pact point  (a  "z  lattice  unit"  is  1.414  times  greater  than  either  the  x  or 
y  lattice  units), 

(2)  x,y,  and  z  velocity  values  in  "z  lattice  units"  per  second, 

(3)  kinetic  energy,  potential  energy,  and  total  energy  in  electron 
volts,  and 

(4)  vreal  elapsed  time  in  seconds. 

Except  for  a  few  6000  cycle  runs  all  of  the  simulations  used  a  2000 
cycle  cut  off. 
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6.    RESULTS  AND  CONCLUSIONS 

The  data  collected  from  the  printout  sheets  discussed  in  section 
five  were  condensed  and  tabulated  as  follows :  (1)  primary  final  energy 
(either  approximately  25ev  or  the  energy  at  the  2  000  or  6000  cycle  cut 
off),   (2)  depth  of  penetration  in  lattice  units,  and  (3)  the  maximum 
spread.    Recall  that  spread  is  the  perpendicular  radial  distance  (in 
lattice  units)  that  a  primary  deviates  from  its  initial  straight  line  tra- 
jectory. 

The  type  of  tracks  reported  by  RO  (those  of  large  penetration  and 
spread)  suggest  the  need  for  a  word  descriptive  of  their  nature.    There- 
fore, in  this  report,  the  term  wander  will  apply  to  all  primaries  that  ex- 
ceed two  lattice  units  of  spread,  but  whose  energy  loss  rate  does  not 
exceed  25ev/A. 

Four  combinations  of  potential  and  energy  were  studied;    Bohr  and 
lOOOev  primaries,  Bohr  and  5  000ev  primaries  ,  Gibson  II  and  lOOOev 
primaries,  and  Gibson  II  and  5  000ev  primaries.    Since  the  data  was 
collected  from  the  impact  point  grids,  it  is  best  displayed  using  a  pic- 
torial representation  of  the  grid.    The  penetration  values  are  shown  at 
their  respective  grid  impact  points  and  can  be  visualized  as  a  penetra- 
tion contour  surface  (figures  9  through  12).    The  pertinent  data  from  the 
spread  contours  was  summarized  and  is  discussed  below. 
A,    Spread  results 

The  spread  values  for  each  impact  triangle  are  a  function  of  primary 
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energy,  potential  and  impact  point.     Separate  spread  contours  are  in- 
cluded in  the  report,  and  are  located  on  the  same  page  as  their  corres- 
ponding penetration  contours. 

The  results  of  the  spread  data  collected  are  summarized  below i 
(1)  Wander,  as  observed  in  the  Bohr  lOOOev  primary  case  by  RO,  was 
found  to  exist  for  several  other  combinations  of  primary  energy  and  poten- 
tial function. 

TABLE  2 


Percent  of  primaries  that  exceeded  2LU  spread 
and  Percent  of  primaries  that  wandered 

POTENTIAL 
Gib  II                                                 Bohr 

Energy 
lOOOev            5000ev 

Energy 
lOOOev              5000 ev 

Percent 
>2LU  Spread 

27 

54 

43 

25 

Percent 
Wander 

0 

11.5 

5 

15 

The  result  of  RO    can  be  generalized °    When  parameters  are  such  that 
penetrations  vary  from  deep  to  short  values  over  the  surface  of  impact 
triangle,  then  some  of  the  primaries  will  wander.    Note  that  the  Gibson 
II  -  lOOOev  series  of  runs  does  not  satisy  this  criteria  (refer  to  table  2), 
since  all  the  primaries  had  short  penetrations. 

From  another  viewpoint,  if  the  strength  of  interaction  between  primary 
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and  crystal  varies  over  the  impact  grid  from  hard  to  soft  collision,  then 
some  primaries  will  wander.    Thus,  the  presence  of  wander  is  not  simply 
a  result  of  the  Bohr  potential,  but  is  a  function  of  the  strength  of  inter- 
action. 

(2)  A  study  of  the  spread  data  (figure  9  through  12)  reveals  that  the 
primaries  that  wander  are  located  in  a  particular  region  of  the  impact 
triangle.    Note  that  the  primaries  impacting  in  zone  A,  figure  14,  have 
large  spread  values  while  those  in  zone  C  have  large  penetrations.  It 
is  not  surprising  that  the  primaries  that  wandered  were  located  at  the 
intersection  of  zones  A  and  C.    The  overlap  of  the  two  zones  varies  some- 
what with  the  parameters  of  the  run  and  cannot  be  firmly  located  .     How- 
ever, the  important  result  is  that  wander  is  a  function  of  impact  position 
and  is  simply  the  transition  region  from  the  short  range    penetration  zone 
to  the  deep  penetration  impact  zone. 

(3)  The  outlined  rectangular  area  in  zone  C  (figure  14)  has  primaries 
that  have  low  spread.    Any  atom  impacting  this  area  will  have  large   pene- 
tration and  will  not  exceed .8Lu  spread.    The  primaries  in  the  remaining 
area  of  the  impact  triangle  have  at  least  lLu  of  spread  except  for  a 

small  area  near  point  11;  the  midpoint  between  the  atom  crystals.    This 
midpoint  is  a  symmetrical  force  region  for  the  primaries  and  little  deflec- 
tion or  spread  takes  place.    Further,  in  the  Bohr  potential  cases,  the  force 
curve  is  near  cut  off  so  that  the  forces  on  the  primary  are  negligible.    Two 
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regions  of  large  spread  occur  between  the  .3  and  .  5Lu    radii  of  the  cor- 
ner atoms.    These  are  the  zones  where  the  primaries  wander. 

(4)  Primaries  2S  5A,  and  8A  in  figure  15  wandered.    As  shown  by 
the  "history  curves"  of  the  runs,  the  rate  at  which  energy  is  lost  is  us- 
ually uniform  until  approximately  1000-1 2 OOev  and  then  the  energy  loss 

is  more  rapid,.    The  general  appearance  of  the  family  of  curves  is  uniform, 
and  one  of  gradual  transition  from  short  to  long  penetration  as  the  impact 
points  move  from  the  corner  atoms  of  the  triangle  to  the  channel  area 
near  the  right  angle.    The  short  penetration  curves  display  less  regularity. 

The  gap  in  the  region  of  180A  to  300A  penetration  is  not  an  amonaly. 
If  primaries  were  concentrated  in  the  region  of  rapid  penetration  change 
(i.e.  the  "hills"  of  the  penetration  contour)  then  this  gap  would  be  fill- 
ed o    Though  the  crystal  effect  on  range  and  spread  seems  logical,  the 
rigid  consistency  and  smooth  transition  is  surprising.    For  example,  if 
more  impact  points  are  studied,  the  occurance  of  occasional  penetrations 

o  o 

in  the  range  of  300A    to  600A    that  do  not  wander  can  be  expected. 

In  brief  summary,  primaries  do  wander  if  the  impact  point  is  located 
in  a  transition  region  from  hard  interaction  to  soft  interaction.    The  "wan- 
der phenomenon"  is  part  of  a  smooth  change  from  one  type  of  interaction 
to  another. 

(5)  The  four  impact  points  of  the  RO  study  that  wandered  were  investi- 
gated (figure  2).    Their  parameters  were  duplicated  (Bohr  truncated  poten- 
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tial,  lOOOev)  to  insure  a  valid  comparison .     In  addition,  the  Gibson  II 
and  Gibson  III  potentials    were  used,    The  selection  of  the  exact  para- 
meters in  the  Bohr  case  provided  a  comparison  between  the  RO  binary 
model  and  the  n-body  model. 

Of  the  four  points,  only  point  C  displayed  wander  and  it  did  so  with 
the  following  combinations  of  paramenters°    (1)  Bohr  and  lOOOev  pri- 
maries,   (2)  Boivr  and  5000ev  primaries,  and  (3)  Gibson  and  5  000ev  pri- 
maries.   As  seen  in  figure  14,  points  A,  B,  and  D  are  in  zone  C  and 
are  typical  of  that  zone  since  they  displayed  deep  penetration  and  less 
than  .8Lu  spread.    The  probable  reason  for  the  disagreement  between 
the  RO  study  and  this  study  is  the  use  of  the  n-body  model.     (Refer  to 
the  next  subsection  for  detailed  comparison  of  the  results  obtained  from 
the  two  models . ) 

At  present,  none  of  these  results  can  be  verified  experimentally.     If 
the  potential  functions  studied  are  not  reasonable,  the  conclusions  dis- 
cussed above  could  be  radically  change.    As  will  be  seen  in  the  next  sec- 
tion, the  potential  functions  selected  (particularly  the  Bohr)  may  be  too 
weak.     If  so,  the  amount  of  "wander"  found  in  future  studies  using  "hard- 
er" potentials  may  decrease  or  disappear. 
B.     Penetration  results 

The  penetration  results  are  shown  in  contour  diagrams  (figure  9  through 
12)  and  the  integral  penetration  curves  of  figure  6.     Penetration  contours 
for  all  runs  are  regular  and  well  defined.    The  following  points  are  reveal- 
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ed  from  the  data° 

(1)  Impact  points  4,  4A,  7  and  7AC  which  are  located  in  the  110 
channel  of  the  crystal  have  maximum  penetration .    The  Gibson  II- 
lOOOev  contour  is  complete  and  reveals  the  smallest  and  lowest 

peak  of  maximum  penetration.    The  Gibson  II-5  000ev  contour  is  virtual- 
ly complete  (four  runs  truncated  at  600Lu)  and  is  shaped  similarly  to 
the  Gibson  II  lOOOev  contours „    However,  the  contour  height  is  much 
greater  and  the  region  of  very  low  penetration  (below  20Lu)  is  much 
smaller. 

Both  Bohr  contours  were  truncated  at  2000  cycles.  In  comparison  to 
the  Gibson  results,  the  areas  of  low  penetration  are  smaller  and  the  rate 
of  change  to  deeper  penetration  more  rapid.    Excessive  computer  run 
time  curtailed  the  construction  of  complete  contour  surfaces  for  the  Bohr 
potentials. 

The  Bohr-lOOOev  contour  and  the  Gibson  II-5  000ev  contour  are  simi- 
lar in  size  and  shape,  which  indicates  that  these  combinations  of  ener- 
gy and  potential  produce  similar  interactions „ 

(2)  As  the  energy  of  the  primary  is  increased  the  depth  of  penetra- 
tion increases;  also,  the  amount  of  penetration  varies  with  the  potential 
function  used. 

The  n-body  integral  penetration  curve  (5  000ev  and  Gibson  II  poten- 
tial has  much  larger  penetrations  than  the  RO  results  (figure  6).     Un- 
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fortunately,  none  of  these  simulation  curves  can  be  easily  compared 
to  the  experimental  data  of  figure  1,  because  of  the  large  difference 
in  the  substances  studied.     Important  conclusions  are  available  from 
the  results,  however. 

The  5  000ev  curves  in  figure  6  of  RO  and  this  study  agree  fairly 

o 

close  in  the  region  of  small  penetration  (up  to  120AK    Thus,  it  is  pro- 
able  that  the  n-body  program  agrees  closely  with  the  binary  model  for 

o 

severe  interactions  (he.   small  impact  parameters  below  .  3Lu  or  .  5A), 
A  thorough  study  by  GH  established  the  close  agreement  between  single 
binary  collisions  and  n-body  collisions  at  short  impact  parameters. 
Their  study  investigated  a  wide  range  of  energies  and  several  poten- 
tials (including  the  Gibson  II) . 

(3)    In  the  region  of  deep  penetration  the  two  models  disagree  . 
Nevertheless,  an  important  fact  was  revealed  when  the  simulation  pene- 
tration data  were  compared  with  the  corresponding  potential  curves. 
Refer  to  figure  13  which  shows  the  Gibson  II  potential  discussed  here, 
the  Bohr  potential,  and  a  new  potential  to  be  explained  in  the  next  sub- 
section.   As  seen  in  the  diagram,  the  three  curves  in  question  use  the 
basic  Gibson  II  potential  function.     However,  RO  curves  B  and  C  are 
truncated  at  much  larger  potential  values  than  the  potential  of  this  study. 
As  seen  from  the  integral  penetration  curve  and  these  potentials,  as 
the  distance  to  the  truncation  point  is  increased,  the  depth  of  penetra- 
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tion  of  the  primaries  increase  .     (This  effect  is  opposite  of  what  might 
be  first  supposed).    The  most  plausible  explanation  probably  lies  in 
the  subtle  effect  of  the  small  potential  curve  slope  (hence  small  force 
values)  that  occured  in  the  region  beyond  the  truncation  of  the  RO 
curves.    These  small  forces  confine  the  primaries  to  the  more  open 
crystal  regions.    As  an  analogy,  the  deeper  penetrations  can  be  com- 
pared to  a  rock  skipping  over  the  surface  of  water „     It  makes  many 
weak  impacts  that  confine  its  movement.    This  is  the  case  of  the  n- 
body  potential.    The  more  radically  truncated  potentials  of  the  RO  study 
probably  cause  abrupt  interaction  when  a  primary  moves  from  a  region 
of  no  potential  to  one  of  large  finite  value.    Thus,  deep  penetration  of 
primaries  still  occurs  in  the  RO  model  but  to  a  much  smaller  extent. 

In  summary,  one  significant  factor  causing  the  deeper  primary  pene- 
trations of  this  study  was  the  use  of  a  potential  that  was  truncated  at 
large  radii.    The  use  of  the  n-body  model  also  contributed  to  large  pene- 
trations and  will  be  discussed  shortly. 

(4)    As  shown  in  figure  13,  the  RO  study  also  used  an  eroded  poten- 
tial  that  decreased  to  zero  at  1.27A.    The  corresponding  integral  pene- 
tration  curve  is  similar  in  appearance  to  the  n-body  penetration  curve  of 
this  study,  but  has  slightly  lower  penetration  values.  Thus,  the  poten- 
tial function  if  eroded  or  truncated  at  large  distances  tends  to  produce 
similar  large  penetration  depths,     In  contrast,  potentials  that  have  signifi- 
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cant  values  at  truncation  reduce  penetration  depth. 

(5)   Another  factor,  which  undoubtedly  caused  the  deeper  penetra- 
tions found  in  this  study,  is  the  n-body  modeL     It's  use  insures  a 
continual  superposition  of  all  significant  crystal  forces;  and  consequent- 
ly, the  net  force  on  the  primary  is  usually  smaller  and  produces  less 
scatter.    Thus,  the  wander  detected  by  RO  for  impact  primaries  A,B,  and 
D,  but  which  was  not  verified  using  the  n-body  model,  may  result  from 
the  insufficient  superposition  of  forces  in  the  binary  approximation  collis- 
ion process.    This  helps  clearify  the  results  in  the  previous  subsection 
concerning  wander. 
C.    Recommendations 

Two  proposed  simulation  studies  are  presented  in  this  subsection. 

Both  are  attempts  to  obtain  verification  with  experimental  data.    Recent 

18 

studies  by  Abrahamson       provide  potential  functions  for  tungsten  and  the 

inert  gases.    This  development,  in  conjunction  with  the  experimental 
work  of  Davies,  provides,  for  the  first  time,  the  possibility  of  positive 
experimental  verification  of  computer  simulation  programs. 

It  is  recommended  that  the    n-body  program  incorporate  the  new  poten- 
tial function  developed  by  Abrahamson,  and  that  it  be  used  to  simulate 

135 
the  Davies  experiment  of  Xe         in  tungsten. 

Although  the  experimental  results  for  copper  primaries  in  copper  are 

not  available,  it  is  felt  that  the  simulation  penetration  values  are  too 
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large  (especially  the  n-body)  and  should  agree  more  closely  with  the 
25  kev  -  Kr        in  Cu  experiments  and  the  Xe         in  Cu  experiments  at 

5kev.    In  this  regard  RO  have  used  "Cu"  primaries  of  triple  mass 

192 

(Cu        )  which  did  reduce  penetration  ranges  (fig.  6),  but  not  suffi- 

84 
ciently  to  agree  with  the  results  of  25kev  Kr       in  Cu. 

In  spite  of  the  probability  that  the  n-body  integral  penetration  curve, 
obtained  in  this  study,  is  further  from  experimental  fact  than  the  RO 
curve,  the  n-body  model  is  an  inherently  more  accurate  approach.    This 
conclusion  is  based  on  previous  discussion  which  reduces  to  the  essen- 
tial fact  that  the  n-body  model  is  more  realistic  and,  hence,  a  better 
simulator  „    Equally  critical,  the  results  of  the  RO  study  were  obtained 
using  severe  truncation  values  (fig  13).    Truncation  of  the  force  func- 
tion at  significant  values  is  unrealistic  and  probably  leads  to  incorrect 
results.    To  obtain  closer  verification  of  experimental  results,  stronger 
potential  functions  may  be  necessary    in  future  n-body  program  studies. 
(It  is  possible  that  experimental  integral  penetration    curves  could  be 
verified  by  most  models  by   sufficient       adjustment  of  the  potential 
functions.    Thus,  experimental  verification  in  a  single  area  will  probab- 
ly not  prove  sufficient,  in  the  long  run,  to  establish  a  particular  poten- 
tial function). 

Preliminary  calculations  indicate  that  the  potential  function  of 
Abrahamson  (Xe  in  W)  is  significantly  stronger  than  the  Gibson  II  poten- 
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tial.     If  so,  it  is  likely  that  this  potential  will  produce    excessive 
duction  in  primary  penetrations  and  its  use  in  a  simmulation  program  will 
not  yield  close  correlation  to  Davies  experimental  work. 

It  is  also  recommended  that  the  verification  of  Davies  experiment 
be  attempted  using  a  "compound"  force  function.    This  function  should 
follow  the  Gibson  II  force  curve  at  small  impact  parameters  where 
good  agreement  is  obtained  at  present.    At  approximately  1.0A    the 
curve  would  change  to  a  different  exponential  form  that    exerts  much 
larger  long  range  forces  on  the  primaries.    This  function  will  reduce 
deep  penetrations  and  agree  more  closely  with  experimental  data.    The 
parameters  of  Davies  Xe  in  W  should  be  duplicated  for  this  study  also. 
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Fig.  1.  INTEGRAL  PENETRATION  CURVES 
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Fig,  2.  ROBINSON  AND  OEN  "...(Oil)  channel  trajectories  onto  (Oil)    i3 
surface  for  f.c.c.  Cu...(lKev  Cu  primaries )... Bohr  potential." 
Points  A,B,C  and  D  penetrated  300A  in  part  of  trajectory  shown. 
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Fig.  9.  CONTOUR  CURVES.  Gibson  II  potential  and  ^OOOev  primaries, 
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Fig,  10.  CONTOUR  CURVES,  Gibson  II  potential  and  lOOOev  primaries, 
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Fig.  12.  CONTOUR  CURVES.  Bohr  potential  and  lOOOev  primaries, 
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APPENDIX  II 
PROGRAM  DISCUSSION 
A.     Program  Explanation 

The  descriptions  given  below  are  to  be  used  with  Appendix  I  (Pro- 
gram Listing),  and  are  specifically  applicable  to  programs  Chan3  ,  Bull- 
ets, or  G  Chan3. 
BOX  1 

Listing  of  dimension,  common,  and  Format  statements.     Provide  graph 
plot  data  needed  such  as  graph  size,  scale,  and  titles.    Read  force  and 
potential  function  parameters.    Read  primary  atom  parameters.    Evaluate 
miscellaneous  constants. 
BOX  2 

Designation  of  nearest  neighbor  lattice  distances.     Construction  of 
microcrystallite.    Three  nested  "DO  loops"  that  extend  to  statement  63 
are  used  to  locate  and  number  the  atoms  in  the  crystal. 
BOX  3 

Initialize  the  velocity  and  the  position  of  the  primary. 
BOX  4 

Store  initial  positions  of  all  atoms.     Initialize  various  terms  (i.e. 
set=0).     Print  variable  quantities  such  as  primary  energy,  direction 
cosines,   mass,  etc.  which  comprise  the  output  heading. 
BOX  5 

Compute  several  constants  given  DT,  CVD,  PTMAS  and  PGMAS. 
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BOX  6 

This  entire  box  is  a  "DO  loop"  on  all  crystal  atoms.     Set  all  forces 
to  Oo     Check  coordinate  distances  from  primary  to  lattice  atoms  and 
neglect  if  greater  than  ROE,     Compute  the  square  of  the  "vector"  dis- 
tances between  the  primary  and  lattice  atoms  and  neglect  if  greater  than 

2 
(ROE)    .     Compute  the  force  that  corresponds  to  the  vector  distance  and 

compare  it  to  FM.     Neglect  if  less  than  FM„     Compute  forces  on  the 

primary  and  lattice  atoms, 

BOX  7 

If  index=l  jump  to  box  8.     Compute  temporary  position  of  the  primary 
and  then  the  rest  of  the  crystal  atoms.    All  atoms  are  advanced  at  their 
initial  velocity  for  DT  seconds,  while  experiencing  forces  found  in  box  6, 
BOX  8 

(Return  to  box  6  and  recompute  the  forces  on  all  atoms.     These  values 
are  added  to  the  initial  forces  yielding  force  terms  twice  the  average. 
Since  index=l  on  reaching  box  7  a  jump  is  made  to  statement  260  in  box  8). 
The  averages  of  the  initial  and  temporary  position  forces  are  used  to  compute 
the  final  positions  of  the  primary  and  lattice  atoms,     (See  page  14  of  text). 
BOX  9 

Perform  tests  to  see  if  program  output  is  necessary  or  if  the  program 
should  rebuild  the  crystal.     If  the  latter  is  true  jump  to  the  appropriate  re- 
generation box. 
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BOX  10 

Calculate  the  distance  from  the  primary  to  each  lattice  atom;  if  less 
than  ROE  calculate  the  potential  energy  of  the  interaction  and  store. 
BOX  11 

Calculate  the  quantities  needed  in  the  output  statement  (i.e.  pri- 
mary position,  kinetic  energy,  total  energy,  and  real    elapsed  time). 
Write  the  output  statement.     Check  to  insure  the  potential  energy  is  less 
than  "Bitty".    Terminate  program  if  the  potential  energy  is  less  than  "Bitty1 
for  three  passes. 
BOX  12 

Check  current  program  cycle.     If  at  "cut  off"  cycle,  terminate  the 
program. 
BOX  13 

This  section  regenerates  the  microcrystallite  in  the+Y  direction. 
BOX  14A 

Regenerate  microcrystallite  in  the  -X  direction. 
BOX  14B 

Regenerate  microcrystallite  in  the  +X  direction. 
BOX  15A 

Regenerate  microcrystallite  in  the  -Z  direction. 
BOX  15B 

Regenerates  microcrystallite  in  the  +Z  direction. 
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BOX  15 

Call  the  graph  plot  subroutine  (called  Draw) .    The  following  parameters 
are  necessary: 


Variables 
II 

GRX 
GRZ 

C 

0 

IABEL 

ITITLE 

EXSCALE 

YSCALE 

0 

0 
2 
2 
9 
15 
1 
ERROR 


Explanation 

Number  of  points  to  be  plotted 

X  coordinate 

Y  coordinate  (simulation  program 
Z  coordinate) 

Number  of  curves  per  graph 

Point  or  curve  plot 

Type  End,  by  the  end  of  curve 

Graph  title 

Graph  scale  in  X  direction 

Graph  scale  in  Y  direction 

Distance  of  x  axis  from  bottom 
of  graph 

Distance  of  y  axis  from  left  edge 

Mode  of  x  axis  (automatic) 

Mode  of  y  axis  (automatic) 

Graph  width 

Graph  heighth 

Place  grid  on  graph 

Indicates  if  previous  plot  was  OK 
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B.     Definition  of  Variable: 

Term 
Bitty 

BE 

BX,  BY,  BZ 

COX,  COY,  COZ 

CVD=1.8075xlO"10 
19 


■27 


CVE= 1.6 0x10 
CVM=1.  672x10" 
DIST 

DT 

DTI 

DTOD 

DRX,  RRY,  DRZ 

DXS,  DYS,  DZS 

DXT,  DYT,  DZT 

EI=25.0 

EXA,  EXB,  EXC,  EX1 

EV 


Usage 

a  small  increment  that  is  compared  to 
primary  potential  energy  (PTE).     Used 
to  terminate  program  when  potential 
energy  goes  below  increment  value. 

primary  energy  (in  joules) 

primary  initial  coordinate  (in  lattice 
units) 

direction  cosines  of  primary 

o 

converts  lattice  units  to  A  (Cu  only) 

converts  ev  to  joules 

converts  AMU  to  kg. 

distance  between  two  atoms  or  the 
square  of  the  distance  (lattice  units) 

basic  time  step  length  (seconds) 

variable  scale  factor  for  DT 

DT/CVD  (a  "units  converter") 

coordinate  d  istance  of  primary  to 
lattice  atoms  (in  lattice  units) 

locates  current  microcrystallite  with 
respect  to  original  origin  (lattice  units) 

test  constants  used  to  rebuild  crystallite 

cut  off  energy  (in  electron  volts) 

input  potential  function  parameters 
(see  Gay16  thesis) 

primary  energy  (in  electron  volts) 
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FA,  FB,  FC 

FM 

FOD 
FORCE 

FXA,  FXC 


GMAS 
GRX,  GRZ 

IDX#  IDY,  IDZ 
IT,  JT,  KT 
K,  T£,  IZ 

dcp,  TfP,  izp 

IYT 

HDTOM 
HDTOMB 
HGMAS 

HTMAS 

M=2 


x,y,  and  z  force  components  be- 
tween two  atoms  (in  newtons) 

a  small  force  increment  used  to  eli- 
minate weak  force  terms 

force/distance  (in  newtons/lattice  units) 

eroded  force  between  two  atoms 
(in  newtons) 

input  force  function  para  meters 
(see  Gay  thesis 

x,y,z  coordinate  force  terms,  respectively 
(in  newtons) 

primary  mass  (in  AMU) 

coordinates  of  primary  that  are  used  by 
graph  plot  routine  (in  lattice  units) 

constants  used  in  lattice  reconstruction 

test  points  in  lattice  generator 

unite  lattice  dimensions  (in  lattice  units) 

specify  volume  of  microcrystallite  (in 
lattice  units) 

used  in  crystal  regeneration 

DT/(2*PTMAS) 

DT/(2*PGMAS) 

V2  primary  mass,   units:  (ev*kg)/ 
(m2/sec2) 

1/2  target  mass,  units;  (ev*kg)/ 
(m2/sec2) 

a  constant 
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MASK 

ND 
NS 
NT 
NTT 

ONE=1.0 

OPNX 

OPNZ 

PGMAS 

PKE 

POT 

PPE 

PTE 

PTMAS 

QDTOM 

QDTOMB 

ROE 

RO 

RX,  RYf  RZ 

RXI,  RYI,  RZI 


used  to  calculate  absolute  value 
(removes  minus  signs) 

print  statement  increment  (in  cycles) 

initial  print  statement  (in  cycles) 

current  program  cycle 

program  cycle  for  final  print  statement 
("cuts  off  program11) 

a  constant 

1.9  SCX 

1.9  SCZ 

primary  mass  (in  kg) 

primary  kinetic  energy  (in  electron  volts) 

potential  between  any  two  atoms  (in  elec- 
tron volts) 

primary  potential  energy  (in  electron 
volts) 

primary  total  energy  (in  electron  volts) 

target  mass  (in  kg) 

DT/ (4*  PTMAS) 

DT/(4*PGMAS) 

nearest  neighbors  separation  (in  lattice 
units) 

=  l/ROE 

atom  coordinate  (in  lattice  units) 

initial  atom  coordinates  (in  lattice  units) 
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RXK,  RYK,  RZK 
SCX,  SCY,  SCZ 

TIME 

TMAS 

TPO=2.1 

TPOX 

TPOZ 

VOL 

VX,  VY,  VZ 

XLL,  ZLL 

ZE  =  0 .  0 


temporary  atom  coordinates  used  in 
force  calculation  (in  lattice  units) 

factors  to  convert  X,Y,Z  to  the  same 
lattice  unit  basis;  (values  depend  on 
type  crystal  used;  i«e0  BCC  or  f.c.c, 
etc.) 

real  elapsed  time  (in  seconds) 

target  mass  in  AMU 

a  constant 

2 . 1  SCX 

2  . 1  SCZ 

velocity  of  primary  (meters/seconds) 

atom  velocity  (in  meters/seconds) 

test  constants  that  determine  if 
regeneration  is  necessary 

a  constant 
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APPENDIX  III 
PROGRAM  OPERATING  PROCEDURES 
The  simulation  model  consists  of  a  single  basic  porgram  that  has 
two  forms: 

(1)  A  fast  running  program  that  can  reconstruct  the  basic  lattice  in 
the  y  direction  only  (the  direction  of  the  incoming  bullet) .    This  pro- 
gram is  limited  to  cases  in  which  the  bullet  does  not  intereact  strong- 
ly and  hence  doesn't  leave  the  basic  lattice  in  any  direction  except 
the  y  direction.    This  program  runs  in  about  one  half  the  time  of  the 
other  program  and  has  three  identical  decks  which  are  titled  as  follows : 
Chan,  Gchan,  and  Bulleto    For  runs  cut  off  as  2000  cycles  the  "Chan 
type  "  program  requires  about  ten  minutes.,    Those  shots  that  strike 

the  impact  triangle    near  the  corner  atoms  often  terminate  in  very  short 
time  since  the  primary  rapidly    drops  below  the  cutoff  energy  of  25ev. 

(2)  The  slower  program  can  contain  any  bullet  trajectory  since  it 
reconstructs  the  lattice  in  any  direction  required.    The  testing  procedure 
to  perform  this  function  increases  the  program  run  time.  This  program 
has  two  identical  decks  that  are  titled;  Chan3  and  BuJet3.     (The  various 
deck  names  serve  only  to  expedite  processing  and  ease  administration 
of  runs.    The  slower  program  (Chan3)  was  the  third  modification  of  the 
basic  n-body  program).    Run  times  are  about  20  minutes  for  2000  cycle 
cases. 

The  items  discussed  in  the  remainder  of  this  section  are  applicable 
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to  all  decks: 

A.  Parameter  cards  internal  to  the  deck,, 

Within  the  program  deck  there  are  three  parameter  cards  that  are  often 
changed;  the  potential  function,  force  function,  and  the  force  truncation 
card.    The  force  and  potential  cards  are  the  third  and  fourth  cards  in  all 
decks.     These  fortran  functions  are  computed  from  a  single  imput  data 
card  (see  the  next  sub  section)  „    Thus  the  force  and  potential  cards  are 
paired  to  an  input  parameter  card.     Specifically,  the  Bohr  potential  and 
force  cards  have  a  single  corresponding  input  cardo    The  Born- Mayer 
force  and  potential  cards  become  one  of  the  three  Gibson  forms  (1,11,  or 
III)  through  the  use  of  one  of  three  corresponding  input  parameter  cards. 
The  third  internal  card  is  the  force  card  (ex.  fm  =  1„0E-11)  which  allows 
the  selection  of  any  desired  truncation  force. 

B.  Input  data  cards . 

These  are  the  basic  data  cards  used  at  the  end  of  the  deck  to  spec- 
ify the  various  parameters  for  a  run  and  to  set  up  the  graph  plot  routine. 
The  following  list  describes  the  cards  in  order  as  they  would  appear  for 
a  single  rune 

(1)  The  first  card  prints  any  desired  statement  at  the  top  of  the 
print  out. 

(2)  This  is  the  force  and  potential  input  parameter  card  and  it  supplies 
the  follwoing  values?  EX1 ,  EXB,  EXC,  FXC,  IH.     The  first  four  quantities 
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are  needed  to  define  the  function  and  the  last  (IH)prints  the  potential 
function  in  the  printout .  The  above  data  cards  are  used  once,  per  run. 
(Chan5  uses  a  new  input  parameter  system). 

(3)    This  card  provides  the  variable  quantities  that  specify  the 
primary  parameters  and  general  program  information.    The  following 
items  appear: 

GMAS  =  Primary   mass  in  AMU 
TMAS  =  lattice  atom  mass  in  AMU 

BX  =  X  initial  coordinate  for  the  primary 
BY  =  Y  initial  coordinate  for  the  primary 
BZ  =  Z  initial  coordiante  for  the  primary 
EV  =  Primary  energy  in  ev 
COX  =  Direction  cosines  of  primary 
COY  =  From  the  Oil  plane  to  trajectory  line 
NTT  =  Cut  off  cycle  (ends  program) 
NS  =  First  print  statement  (in  cycles) 
ND  =  Print  increment  (in,  cycles) 
Recall  that  one  cycle  =  .  1  time  for  the  primary  to  travel  one  lattice  unit 
(1.807A). 

The  next  six  data  cards  are  used  to  title  and  identify  the  graph  out- 
put for  the  run. 
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(4)  Potential  is  designated  (ex.  GIB.  for  Gibson) 

(5)  Energy  (ex.  5  000ev) 

(6)  X  coordinate  (ex.  x  =  2.  00) 

(7)  Z  coordinate  (ex.  z  =  3.00) 

(8)  Cosine  X  (ex.  cox  =  .00) 

(9)  Cosine  Y  (ex.  coy  =  1.00) 

This  completes  the  data  cards  needed  for  single  run. 

To  provide  additional  runs  cards  three  through  nine  are  repeated  for 
each  run.    An  "end  card"  is  placed  last  on  the  deck  to  terminate  the  pro- 
gram.   The  program  is  designed  to  run  on  the  regular  monitor  operation 
service,  but  for  runs  longer  than  15  minutes  (the  operation  service  limit) 
the  program  must  be  run  individually.    Three  tapes  are  required,  in  add- 
ition to  the  Fortran  60  compiler  tape,  to  run  the  program  on  the  1604 
computer.    These  tapes  are  the  input,  output,  and  graph  plot  tapes. 

(To  perform  runs  without  the  graph  plotter  pull  the  "call  graph  card", 
the  title  cards  (4  thru  9),  and  the  card  that  "reads  in"  the  title  card. 


72 


